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In connection with the design of transistor circuits, for example, it is 
frequently necessary to obtain a numerical solution of a system of non- 
linear ordinary differential equations. In some cases, these equations 
possess a property that leads to intolerable computational requirements 
relative to the use of standard predictor-corrector techniques or general 
linear multipoint formulas of open type. 

Here we describe an alternative approach which has been used to solve 
some practical problems by permitting dramatic step-size increases (for 
example, a factor of 10*). The approach is developed in a way which 
provides some detailed understanding of why it is useful. 

I. INTRODUCTION 

In connection with the design of transistor circuits, for example, it 
is often necessary to obtain a numerical solution of a system of non- 
linear differential equations 

x + f(x, = 0, t = 0, [:r(0) = x ] (1) 

in which x and f(x, •) are N- vector- valued functions of t. The sim- 
plest numerical-integration formula which can be in principle used 
for this purpose is Euler's formula: 

?/ n+ i = y n + hy' n , n ^ (2) 

in which h, a positive number, is the step size; y = x ; 

Vn = —f(yn , nh) for n ^ 0; 

and y n is of course the approximation to x(nh) f or n =• 1. 

It is frequently the case that /(re, •) possess a property that leads 
to computational requirements consistent with the use of (2) that are 
intolerable. To see clearly how this situation can arise suppose that 
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the solution of (1) is desired over some finite interval [0, r], and con- 
sider the very special case in which f(x, t) = Ax with A an N X N 
matrix possessing distinct eigenvalues {a ( } all of which have positive 
real parts. Then using the fact there exists a nonsingular transforma- 
tion T such that 

A = TDT~\ D = diag (a t , a 2 , • • • , a N ) (3) 

we have 

y n+1 = T(1 N - hD^y* , n ^ 0, [y = x ] (4) 

in which 1 N is the identity matrix of order N. From (4) 

y k = T(1 N - hD) k T- x x Q , h^O. (5) 

Since 

x(kh) = Te- Dkh T- l x k ^ (6) 

it is evident that the numerical solution is "acceptable" if h is so small 
that (1 — hai) k is an "acceptable" approximation to e~ aikh for all i 
and all values of k for which ^ kh ^ t. On the other hand if for 
some value of i 

| 1 - ha { | = 1, or | 1 - ha t \ > 1 

then for at least one initial condition vector x , { || y k \\ } o (II • II denotes 
the usual Euclidian norm) does not approach zero as k —» °o or is 
unbounded, respectively [that is, (2) is numerically unstable]. Therefore 
if the sequence \y k ) defined by (4) is to be a good approximation to 
the samples of the solution of (1) with j{x, t) = Ax, it is certainly 
necessary that 

| 1 - ha t | < 1 for all i. (7) 

Moreover, in order to fully determine the character of the solution of 
the differential equation, it is reasonable to assume that t, the length 
of the interval over which the solution is desired, is proportional (by 
some factor c such as 3 or 10) to the reciprocal of min t Re (a<) (that 
is, proportional to the largest time constant of the system) . Thus in 
addition to (7) we have 

r =c[minRe(a,-)]" 1 . (8) 

A lower bound on the number of evaluations of (2) necessary to 
compute the solution is r/h where h satisfies (7) . If all of the a t are 
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real, the smallest lower bound is simply 

max (a,) 



\c '■ , r - (9) 

mm (a,) 

It is a simple matter to give examples of, for instance, positive-element 
linear RC networks governed by a state equation of the form x + Ax = 
for which the bound (9) can be made arbitrarily large by choosing the 
value of one capacitor to be arbitrarily small. Thus, from the practical 
viewpoint, computation based on (2) can be impossible as a result of 
the presence of parasitic circuit elements that have no really signifi- 
cant effect on the circuit performance! It is not surprising therefore 
that a more complex and pressing problem of the same type arises in 
connection with the numerical solution of the nonlinear differential 
equations of transistor circuits, as a result of, for example, the para- 
sitic capacitors associated with the models of transistors. For many 
practical circuits of this type, computation time estimates, based upon 
use of (2) and a modern high-speed computer, are about 1000 hours. 

The well-known basic problem described above arises not only in 
connection of the use of (2), but (as can easily be shown) is en- 
countered also in attempts to use more general integration formulas 
of open type 1 - - 

y v 

y« + i = Z akVn-i + h £ b&. t , (10) 

*-0 t-0 

or predictor-corrector techniques 1, - such as 
2/n +1 (p) - y n -r + 2hy' n 

</ n+1 (c) = !/. + J%i + !/i« w ). 
The fundamental difficulty associated with the integration of "stiff 
equations" results from the restrictions that must be imposed on h in 
order to insure numerical stability. 

The purpose of this paper is to consider the properties of alterna- 
tive numerical methods for obtaining solutions of equations of the 
form (1). Our principal objective is to present some analytical results 
that shed some light on the properties of a class of numerical-integra- 
tion techniques that have been used to solve practical transistor circuit 
problems by permitting dramatic step-size increases (for example, a 
factor of 10*) relative to the methods defined by (10) and (11). 

More explicitly, attention is focused on "l&rge-h algorithms" based 
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on, or derived from, the standard formula of closed type 

Z/» +1 - V* + %n + i ( 12 ) 

which is a special case of the general multipoint formula 

V V 

y« + i = Z a kyn-k + h 2J bkVn-k ( 13 ) 

h=0 k--l 

with 6_! 5^ 0. There is an extensive body of information concerning 

(12) in the numerical-analysis literature only for the case in which 
h is "sufficiently small." 

II. INTEGRATION FORMULA 

If we use the numerical-integration formula 

tf.+i = Vn + hy' n+l (12) 

in an attempt to compute the solution of (1), then y n+1 is denned 
implicitly in terms of y„ through 

y n+ i + hf[y n+1 , (n + l)h] = y n , n ^ 0, [y = x ]. (13) 
For the special case considered in Section I, in which f(x, t) = Ax 
and A = TDT~ l , we have 

y n+1 = T(1 N + H>)- 1 r- , y. , n ^ (14) 

and to the extent that (1 + Aa,-)" 1 is a good approximation to e"° ,A , 

(13) generates an acceptable numerical solution of the differential 
equation. More explicitly (13) generates the exact solution of the 
differential equation 

x + Ax = 0, t ^ 0, [s(0) = sj (15) 

in which A = TDT' 1 and e~ bh = (1 N + hD)~\ 

Let us suppose now that all of the a, are real and that ha t is very 
small relative to unity for i belonging to a proper subset S of 31 = 
{1, 2, • ■ • , N\, and_that Aa, is very large relative to unity for i belong 
to the complement S of S with respect to 91. Then for all i t S, a,- , the 
tth element of D is very nearly a,- , while for all * c S, a, < a, and a,- 
is very much larger than all of the a, for which % e S. 

In other words, roughly speaking, (13) generates a solution to a 
differential equation governing a system similar to that governed by 
x + Ax = 0; the former system has virtually the same low-frequency 
performance and less pronounced high-frequency performance. To 
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look at the situation in still another way, in using (13) we are able 
to (i) break away from an extremely restrictive requirement on h for 
numerical stability, such as (7), and (ii) trade step-size tor accuracy 
of high-frequency solution components. 

The simple heuristic argument given above suggests that the use 
of (12) can lead to a considerable increase in permissible step sizes 
for a class of nonlinear transistor circuit problems in which typically 
the Jacobian matrix df(x, t)/dx of f(x, t) along the solution of (1) 
possesses only real eigenvalues which are widely separated. This argu- 
ment is supported by a proposition, proved in Section IV, which is 
concerned with the case in which there exists a constant w > such 
that (with (• , •) denoting the usual inner product) 

(V, KV, nh) - /(0, nh)) ^m\\y || 2 (16) 

for all w ^ and all y. If this condition is satisfied for all h > 0, which 
for the scalar case is true if 

« a m 

dy 

for all t and all y, if |j /(0, t) || -*■<) as * -» « or if || /(0, t) \\ is uniformly 
bounded on [0, <»), then (as can easily be shown) \\x{t) \\ — > as t — > » 
or || x (t) || is uniformly bounded on [0, co), respectively. The Proposition 
asserts that if (16) is met and y n+1 is defined f or n ^ by (13) ,then 

111/- II ^ (1 + mhy \\x || + £(l + m/i)- ( * +n \\hj[0, (n-k)h] \\ 

k=0 

for all n ^ 1, which implies that (13) is numerically stable for &ll h 
in the sense that for all h, \\ /(0, nh) \\ -^ as n -> <x> implies that y n -*■ 
as n -» oo and {\\ /(0, nh) \\\™ bounded implies that {y n } n is bounded. 

Although the result stated above does not provide quantitative in- 
formation concerning the errors incurred in using (13), it does show 
under a reasonable assumption concerning f(x,t) that unlike all for- 
mulas (10) of open type and unlike predictor-corrector methods such 
as (11), (13) defines for any step size a sequence {y„} which is con- 
sistent with either or both of two possible basic properties of the true 
solution. 

The discussion above does not take into account the fact that at 
each step errors are inevitably introduced in solving the equation 

Z/ n+1 + hj[y n+i , (n + l)h] = y n (17) ; 
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for y n +i. Consider the result of using the iteration scheme 

ifij" - ft ~ hf[y ( n k + \ , (n + l)h], y^i = ft 

which is the usual method described *• - in connection with the theory 
of integration formulas of closed type. For the linear case [that is, 
for /(re, t) = Ax], 



= TE (-hDyTSj n . (18) 



Therefore, if y x denotes the approximation to y x computed from y 
after fa iterations, and if y 2 denotes the approximation to y 2 computed 
from y x after k 2 iterations and so forth, then 

y K = TQ kK @ kK _, ■ • • e t fi tl T l y 
in which 

e*, = diag (i: (-hay, • • • , E (-ha N y)- 

Since (assuming now that all of the a, are real) 



i, (-**)* 



> 1 (19) 



provided that ha { > 2 and k, ^ 1, if fta,- > 2 for some *, then J| y k || -> » 
as A; — » oo for some initial condition y ; independent of the sequence 
fa , k 2 , • • • . Therefore the usual iteration method will reintroduce 
the numerical instability for insufficiently small h which it is our objec- 
tive to avoid.* 

Let us consider now a different and more general approach of solving 
(17) for y n+l . Assume that there exists a positive constant I such that 
1(y, nh) satisfies the Lipshitz condition 

II /(ft , nh) - /(ft , nh) || £ 1 1| ft - y 2 || 

for all n ^ and all y t and y 2 . Suppose also that the smallest eigen- 
value of the symmetric part of the Jacobian matrix df(y, nh)/dy of 

♦Similar instability results for the nonlinear case can be proved. But since 
this is hardly surprising, we shall not consider the matter further. 
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f(y, nh) is bounded from below by m, a positive constant, for all y 
and all n. 
Ideally, we would like to determine the sequence {y„\ " denned by 

Vn+\ + hj[y n+1 , (n + l)h] = y n , n ^ 0. 

Suppose that we determine instead a sequence {y n }t such that 
Vo = y 

\\ Vn - v* \\ ^ e 

and 

2/»* + , + hf[y*^ , (n + l)h] = y n 

for n > in which e is an arbitrary positive constant independent of 
n. In other words, suppose that at each step the local error in solving 
for y n+1 is at most e. Then, according to Theorem 1 (Section IV) 

II ft ~ Vn || ^ «(1 + /(/)(1 + /im)-' 2 (1 + Am)~* 

A = 

for all n > 1, which of course implies the uniform bound 

li Vn - y n II ^ 6(1 + hl)(hm,y l , n fc 1. (20) 

Our assumption concerning df(y,nk)/dy implies that the condition 
(y,Ky,nh) -f(0,nh)) ^ m\\y\\ 2 

of the Proposition is met. Thus it follows from the Proposition and (20) 
that if the local error in solving for y n+1 is held to within e at each step, 
then the algorithm is numerically stable for all h in the sense that 
for all h (i) {|| /(0, nh) ||}* bounded implies that {y„}~ is bounded, and 
(**) II /(0, nA || -> as n -> co implies that for any 8 > there exists 
an n such that || y n \\ ^ e (l + hlKhm)' 1 + 5 for all n ^ n . 

The combination of this stability result and the heuristic argument 
of Section I strongly suggests that the following approach should per- 
mit the use of considerably increased step sizes with acceptable accu- 
racy, for many of the "widely-separated eigenvalue" problems de- 
scribed earlier. Referring to (17), solve for y n+1 at each step using, 
say, the Newton-Raphson technique;* iterate until some norm of 

* After the work reported here had been completed, A. N. Willson, Jr. brought 
to our attention a preprint of a paper by R. Willoughby and several of his col- 
leagues at IBM, in which an approach of this type is suggested. The preprint 
does not contain the principal results of this paper, the material of Section IV 
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the difference between the last two iterates is not greater than some 
small prescribed constant. 

In particular, notice that for f(x, t) = Ax, this approach, using the 
Newton-Raphson iteration procedure, reduces to the use of the for- 
mula y„ +1 = (1 N + hA)~ 1 y n (that is, to equation 14). 

The technique described above has provided a significant reduction 
in total computation time for several types of practical problems. It 
was used, for example, to solve the system of differential equations 
governing the circuit of Fig. 1, an oscillator designed to supply a 1 
kc signal. The 16 G Western Electric 100 Mc. silicon transistor of 
Fig. 1 was represented by a charge-control model (see Section 6.2, 
pp. 556-557 of Koehler 8 ) using two nonlinear charge-controlled voltage 
sources, with the result that the system of equations for the circuit 
is of order 5. 

Motivated by the fact that the local-truncation error for formula (12) 
is ih a x(£) for some £ t[nh, (n + l)h], the following method was used 
(for this problem as well as for others) to control the step size. Let e 
denote the largest of the magnitudes of the elements of the vector 
of second differences associated with the most recently computed point. 
If e e [iS, e] (for this problem e was taken to be 10~ 4 ), then the point 
is accepted; if e > e, then the point is rejected and the calculation 
is repeated with h replaced with %h. lie < \e, then the point is accepted 
and h is replaced by 2h in the computation of the next point. Average 
step-size increases of about 10* (relative to, for example, the use of a 
forth-order predictor-corrector method) were obtained for this problem 
(see Fig. 2). 
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Fig. 1 — One-kilocycle oscillator using a 16G "100 megacycle" transistor. 
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Fig. 2 — Comparison of computed and experimental response of the oscillator 
shown in Fig. 1. 



III. AN EXPLICIT INTEGRATION FORMULA 

Of particular interest in connection with the approach described 
above is the numerical-integration formula 

fr+i - Vn - {l.v + hf[y n , (n + l)h]\- l hf[y n , (n + l)h], 

n ^ 0, [y = x ] (21) 

which is obtained from 

K n+1 + hf[Y n+1 , (n + l)h] = Y n (22) 

by replacing Y n by y n and using as the approximation y n+1 to 
Y n +i the result obtained by using one step of the Newton -Raphson 
iteration scheme with y n the initial point. That is, with 



(23) 



Q(z) = z + hj[z, (n + l)h] - y, , 

</ n+1 - Vn -W(jr) UJ-^C^U. . 

In spite of its relative simplicity, it has been found that formula 
(21) is useful for solving problems of the type that we have been 
considering. For the problem of Fig. 1, it has led to an average step 
size increase of about 10 3 . 
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In view of the simplicity of formula (21), and especially in view of 
the fact that y n+i is defined explicitly in terms of y n , it deserves 
special consideration. 

Theorem 2 (Section IV) asserts that for any h > there exist 
positive constants fc, and k 2 such that k x < 1 and 



|| y n || ^ K || i/o || + hk 2 £ fc} || /[0, (n - fc)ft] || 



(24) 



m 



for n ^ 1, provided that the Jacobian matrix df(y, nh)/dy satisfies 
certain conditions. For the scalar case, these conditions reduce to: 
(i) there exist positive constants k and m such that 

dj(v, nh) , 
dy ~ 

for all y and all n > 1 

for all ?/, n ^ 1 and a e [0, 1]. 

Clearly, under these conditions, y n -> as n -> oo if /(0, nft) -> as 

n -> oo and {?/ n } is bounded if {||/(0, nh)\\} is bounded. 

The function f{y, nh) of Fig. 3 is one for which conditions (i) and 
(ii) are clearly met. If condition (n) is not met, then (24) need not 
follow. To show this, consider, for example, the function of Fig. 4 



--- SLOPE = 1 




Fig. 3 — Definition of f(y, nh) for all n. 
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which meets condition (i), but not condition (u). We have from (21) : 

h = 1 and ?y = 1 imply that y l = —\ 
and 

h = 1 and v/i = — 1 imply that ?/ 2 = 1 
from which it is clear that for this function y n = (-l) n if h = 1 and 
Vo = 1, which of course implies [here /(0, nfr) = for all n] that (24) 
is not satisfied. Thus we see that if condition (ii) is not met, then (24) 
need not follow. 
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Fig. 4 — Alternate definition of f(y, nh) for all n. 
IV. PROPOSITION AND THEOREMS* 

Proposition: If {y n } satisfies 

y»+i + hf[y n+1 , (n + l)h] = y n , n ^ 
and if there exists an m > such that 

(y, 1(y, nh) - 1(0, nh)) £m\\y \\ 2 , n^O 
for all real y, then 

II y n II ^ (1 + mhy n || y n || + E (1 + mh)- k \\ hf[0, (n - k)h] \\ 

t = 

for n > 1. 



♦Throughout this section, ||-|| denotes the usual Euclidean Norm and (•,•> 
denotes the corresponding usual scalar product. 
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Proof: Clearly, 

(y n+ i , Vn - hj[0, (n + l)fe]> 

= <2/„+i , ft»i + hj[y n+l , (n + l)fc] - hf[0, (n + 1)*]> 

£ (1 + mh) \\tf»*i ll 2 > 
and, by the Schwarz inequality, 

(y. + , , V* ~ hf[0, (n + 1)*]> ^ || 2/„ +1 || -II V. II 

+ ||y. + i||-HMO,(n + l)A]||. 

Thus 

II tf„« || Si (1 + mft)" 1 || * || + (1 + m^)" 1 || */[0, (n + 1)*] II 

from which we have 

1 1 y n || £ (1 + mft)"" || y II + Z (1 + w/i)- (t+,) || fc/[0, (n - fc)fc] II 

t = 

for n ^ 1, which completes the proof. 

Definition: Let X(y, nA) denote the smallest eigenvalue of the symmetric 
part of df(y, nh)/dy. 

Theorem 1: Suppose that there exists a constant m such that \(y, nh) ^ 
m > Ofor alln ^ and all y, and that there exists a constant I such that 

II f(Vi , nh) - f(y 2 , nh) || £ Z |l V* " V* II 
for alln ^ awd a/Z t/x and y 2 . If [y n ] satisfies 

if, totiA e a positive constant, \y n ] satisfies 

\\ Vn — y* || ^ e /or 71 ^ mU 
2/* +1 + hf(y* +1 , (n + l)h) = £ n 
Z/ien 

llfc -null £Q + A»)-||fc-* II 

+ (1 + M _1 (l + hl ) * S (1 + ^W 1 )"* 
forn^ 1 . 
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fc+i + fc/BUi + (?/ n * +1 - p n *i), (n + l)h] = y n + (y n+l - y* +1 ) 
and 

y n+i + /i/[2/ n+1 + (y* +1 - $ n+l ), (n + l)h] 

= Vn + hf[y n+1 + fo* +I - £ n+1 ), ( n + i)ft] _ ^[y B+a j ( n + 1 ) fc ] i 
Therefore 

jWi - y n+1 + A/rjf, +1 + (y* +i - y n+1 ), ( n + l)h] 
~ hf[y n+l + (y* +l - y n+1 ), ( n + l)h] 
= y n - y n + (y n+1 - y* +1 ) + hf(y n+1 , (w + l)h) 
- M[y n+1 + (y* +l - y n+1 ), ( n + l)h]. (25) 

With /; the symmetric part of dfty, nh)/dy, the inner-product of 
(&+i - 2/ n+ with the left side of (25) is 

II 0.« - ?/ n+I || 9 + h(y n+l - ? y n+1 , jf /.'{a[£ n+1 + fo* +1 - £ n+1 )] 

+ (1 - «)[?/ n+ . + 0/* +l - y. +1 )], (n + \)h\ da(y n+1 - y n+l )S , (26) 
since 



V=a\ ] + (l-a)[ J 



da(y n+l - y n+i ). 



1\$« + i + (y* +1 - y„ +i ), (n + l)h] - j[y n+1 + (y* +1 - £ n+1 ), (n + 1)A] 

= f df[y, (n + Dh] 
Jo dy 

Expression (26) is hounded from below by 

(1 + hm) \\ y n+1 - y n+1 || 2 . 

By the Schwarz inequality, the inner-product of (y n+l — y n+1 ) with 
the right side of (25) is bounded from above by 

II &+i - 2/ n+ i ||- 1| ft, - y n || 

+ II&+1 - y.+, ||-|| y n+i - y* +1 || + \\y n+1 - y n+1 || 

• II M[y n+1 , (n + l)h] - hf[y n+1 + (y* +1 - y n+l ), ( n + \)h] \\, 
which is bounded from above bv 



^ n+ i - y n+l U-\\ y n - y n \\ + ||ft +1 - y n+1 \\ (« + « 6 ). 
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Thus, 

|| fL»i - </n + i II ^ (1 + *»)"' II fc - 2/n || + (1 + M _1 (l + hl)e, 
from which it follows that 

|| % - Vn || S (1 + M"" || ft, - 2/o || 

n-l 

+ (1 + km)' 1 (I + M) « E (1+ &'»)" 

t-o 

for all n ^ 1. 

Theorem 2: If {y„} satisfies 

y n+ i = y» - [U + *fl* , (» + 1)*]}"%- . (» + IW 

/or n ^ 0, t/ 

(t) ftare ext'ste a constant k < °o swcA $af 



a/(y, gft) 
ay 



< fc 



/or alln ^ 1 and aZZ y 

(ft) tfiere existe a constant m > swcft tfczt X(y, nfc) ^ m for alln ^ 1 
and aZJ y 

(m) irift i^ = hf[y, (n + 1)A] and F a A ftf [«j/, ( n + 1)A], ifce sym- 
metric part of [(2F — F a )Fl\ is* nonnegalive definite for all y, 
all n, and all a c [0, 1], , 

then there exist positive constants fci and k 2 such that &, < 1 and 

|| y„ || ^ fcT II yo II + M, g # II /[0, (n - fc)fe] II for all » £ 1. 
Proof: We have 

y» +1 - y n - { 1* + hf'[y n , (n + l)fc] \- x \hf[y m , (n + 1)A] - MO, (n + 1)*] } 

- {U + hf[y n , (n + l)h]\- l hf[0, (n + 1)A]; 



hence 

llfc+tll* 



•'0 



(/a 



•Il?/n 



+ ||(li, + F)" l ||-||*/[0,(n + l)fc]|| (27) 



* The superscript t denotes matrix transposition. 
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with the understanding that F and F a are evaluated at y = y„, since 

hj[y n , (n + l)h] - hf[0, (n + l)h] = f hj'[ay n , (n + l)h] da y n . 

Jo 

We now prove that there exists fcj c (0, 1) such that 



1* - (I* + F)- 1 f F a da 

•'0 



£*I 



for all n and all y n . 
From condition (Hi) , with 7 an arbitrary iV- vector, 



((2F' - F'„)7, F' a 7) ^ 



or 



(2F'V, F' a V) - (F' a V, F' a V) ^ 
which implies that 

i« 17 112 o/i?'t/ w« tta I II K?«Tr IIS ^ || J7<Tr 112 



F' a V IT - 2(F'7, KV) + II F'V || 2 g || F'V 



or 



|| (f - F'„)7 || 2 g ||F'7 || 2 . 

In view of conditions (i) and (it), it is evident that there exists a 
£ e (0, 1) such that 

-2(FIV, V) £ - (1 - || 7 || 2 - 2(1 - 00*7, V)-(l- || F'7 || 2 
for all «, n, i/ n , and 7. Therefore 

|| (F' - F' a )V || 2 - \\F'V || 2 - 2(FIV, V) 

* -d - & \\ V || 2 - 2(1 - $)(F'V, V) - (1 - || F'7 || 2 
which is the same as 

|| 7 || 2 + || (F' - FDV \\ 2 + 2((F* - F' a )V, V) 

*k\\ 7 || 2 + 2^7, V) + S\\F t V\\* 



or 



|| (1* + f « - Fl)V || 2 =s £ || (1„ + F')V || 2 . 
With U = (lx + F>) 7, we have 

|| (1„ + F l - OOir + J")"'? II 2 £ € || Z7 || 2 . (28) 
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Since (28) is satisfied for all U, 

|| (1„ + F - FDiU + FT 1 II £h 
with k x = U) *. However, 

ll (i* + r- wo, + o _1 II = ll (i* + *)"'(!» + *■ - Will 

and 

1„ - (1„ + J)" 1 f F da 

Jo 

£ f II Uat + F)~\U + F - F a ) || da £ h t . 

Jo 

Consider now 1 1 (In+F^W. Since for any V 

II (1* + F)V || 2 = || V \\ 2 + 2(FV, V) + || FV || 2 £ (1 + 2M || F || 2 , 
it follows at once that 

II (i* + f)- 1 || s (i + atat)" 1 . 

Thus with fc 2 = (l+2ftra)-i 

|| *♦, || ^ h || 2/n || + k 2 || MO, (n + l)h] || 
from which we obtain the bound on \\y„\\ stated in the theorem. 

V. FINAL REMARKS 

The algorithm described in this paper is a marriage of two stand- 
ard techniques, the use of a well-known closed-type numerical-inte- 
gration formula and the Newton-Raphson iteration procedure. It is 
clear that the approach is of use in connection with a certain class of 
practical problems, and, what is of at least as much importance, we 
have some detailed understanding of why the algorithm is useful. 

It is also clear that some natural generalizations and extensions 
of the approach, such as using different closed-type formulas* or 
different methods of solving systems of nonlinear equations, will lead 
to more efficient techniques. Finally, since there are several alternate 
approaches available which are also of use in certain cases (see Pope, 

* For example, for the trapezoidal rule y„+i = y n + \h{y n ' + y n+ i') and for/(x, t) = 
Ax, we have y n+ i = TaT^yn, in which S = diag [(2 - hai)(2 + ton)" 1 , . . . , (2 - 
ha N )(2 + /wAr) -1 ](r and the a { are defined in Section I). In view of the relation 
between the local-truncation errors of the trapezoidal rule and formula (12), this 
suggests that for nonlinear problems the trapezoidal rule should permit larger step 
sizes for the same accuracy when the "fast components" of the solution have decayed 
to a very low level. 
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for example) 4 much work directed toward the comparison of avail- 
able methods is needed. 
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